/*
 *  worker.cpp
 *  phaseit3d
 *
 *  Created by Anton Barty on 27/07/11.
 *  Copyright 2011 Anton Barty. All rights reserved.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "phaseit3d.h"

Worker::Worker(cData3d *diffractionVolume, cData3d *supportVolume, cData3d *guessVolume, cSettings *settings){
	_diffractionVolume = diffractionVolume;
	_supportVolume = supportVolume;
	_guessVolume = guessVolume;
	_settings = settings;
}

void Worker::phaseIt(){

	///	A few arrays for keeping track of things
	long	nx = _diffractionVolume->nx;
	long	nn = _diffractionVolume->nn;
	char	*dvMissingData = (char *) calloc(_diffractionVolume->nn, sizeof(char));;
	char	*support = (char *) calloc(_diffractionVolume->nn, sizeof(char));
	char	filename[1024];


	if(dvMissingData==NULL || support==NULL) {
		printf("Error allocating temporary arrays\n");
		exit(1);
	}

	detectMissingData(dvMissingData);
	detectSupport(support);

	cFFTWdata model = createModelWithRandomPhases(dvMissingData);

	// Start off in Real space
	model.fft(FFTW_FORWARD);

	//	Memory of last guess (for HIO, etc)
	cFFTWdata	lastModel;
	lastModel.create(nx);
	lastModel.createFFTWplan(_settings->nthreads);

	for(long p=0; p<nn; p++) {
		lastModel.data[p][0] = model.data[p][0];
		lastModel.data[p][1] = model.data[p][0];
	}


	/*
	 *	Temporary arrays needed for doing stuff later on
	 */
	cData3d		tempVolume;
	cFFTWdata	tempModel;
	cFFTWdata	shrinkwrapSum;
	cFFTWdata	runningSum;
	tempModel.create(nx);
	tempModel.createFFTWplan(_settings->nthreads);
	shrinkwrapSum.create(nx);
	runningSum.create(nx);
	tempVolume.create(nx);

	for(long p=0; p<nn; p++) {
		shrinkwrapSum.data[p][0] = 0;
		shrinkwrapSum.data[p][1] = 0;
		runningSum.data[p][0] = 0;
		runningSum.data[p][1] = 0;
	}

	/*
	 *	Actual iteration loop
	 */
	time_t			startt, endt;
	clock_t			start, end;
	float			dt;

	time(&startt);
	for(long i=0; i < _settings->niterations; i++) {

		char	processingString[100];

		start = clock();

		/// Start off in real space

		/*
		 *	Calculate real space error
		 *  Error is the amplitude outside of support bounds, divided by total amplitude of solution.
		 *	Effectively % of amplitude outside of the support after Pmod but before applying real space constraints
		 *	Calculated by reference to intensity of solution (saves having to take the square root)
		 */
		float	temp1=0;
		float	temp2=0;
		float	error=0;
		for(long p=0; p<nn; p++) {
			temp1 = CPLX_SQR(model.data, p);
			temp2 += temp1;
			if(support[p] == 0)
				error += temp1;
		}
		if(error != 0.0)
			error = error/temp2;
		else
			error=0.0;


		/*
		 *	Periodically do error reduction
		 */
		int	thisMethod = _settings->prAlgorithm;

		if(_settings->prErrorReductionFreq != 0 ){
            if ( (i-1)%_settings->prErrorReductionFreq == 0 || (i-2)%_settings->prErrorReductionFreq == 0)
                thisMethod = 1;
        }



		/*
		 *	Apply real space constraints
		 *	Object must exist within known support constraint in real space
		 *		1 - simple error reduction
		 *		2 - RASR error reduction
		 *		3 - Fienup hypbrid I/O algorithm
		 */

		switch(thisMethod) {

			// Error reduction
			case 1:
				strcpy(processingString, "err");
				for(long p=0; p<nn; p++) {
					if(support[p] == 0){
						model.data[p][0] = 0.0;
						model.data[p][1] = 0.0;
					}
				}
				break;

			// RAAR
			case 2:
				strcpy(processingString, "raar");
				for(long p=0; p<nn; p++) {
					if(support[p] == 0){
						model.data[p][0] = (1.0 - 2.0*_settings->prBeta)*model.data[p][0] + _settings->prBeta*lastModel.data[p][0];
						model.data[p][1] = (1.0 - 2.0*_settings->prBeta)*model.data[p][1] + _settings->prBeta*lastModel.data[p][1];
					}
				}
				break;

			// HIO
			case 3:
				strcpy(processingString, "hio");
				for(long p=0; p<nn; p++) {
					if(support[p] == 0){
						model.data[p][0] = lastModel.data[p][0] - _settings->prBeta*model.data[p][0];
						model.data[p][1] = lastModel.data[p][1] - _settings->prBeta*model.data[p][1];
					}
				}
				break;
		}


		// Remember current real space array for next iteration
		for(long p=0; p<nn; p++) {
			lastModel.data[p][0] = model.data[p][0];
			lastModel.data[p][1] = model.data[p][1];
		}



		/*
		 *	Operations performed in Fourier space
		 */

		// Go to Fourier space
		model.fft(FFTW_BACKWARD);
		for(long p=0; p<nn; p++) {
			model.data[p][0] /= nn;
			model.data[p][1] /= nn;
		}


		/*
		 *	Fourier space error metric (R-factor)
		 *	Convention is to calculate R from complex amplitudes (Fobs)
		 *	Fobs is defined as amplitude of current guess (and not the intensity).
		 *	Can only be calculated at points where the input data is defined (data bit set).
		 */
		temp1 = 0;
		temp2 = 0;
		float	derror=0;
		float	total=0;
		for(long p=0; p<nn; p++) {
			if(!dvMissingData[p]) {
				temp1 = CPLX_ABS(model.data, p);
				temp2 = _diffractionVolume->data[p];
				derror += fabs(temp2-temp1);
				total += temp2;
			}
		}
		derror /= total;


		/*
		 *	Apply Fourier space constraint
		 */
		for(long p=0; p<nn; p++) {
			if(!dvMissingData[p]) {
				temp1 = CPLX_ABS(model.data, p);
				if (temp1 < _settings->prEpsilon*_diffractionVolume->data[p])
					temp1 = _settings->prEpsilon*_diffractionVolume->data[p];
				model.data[p][0] *= (_diffractionVolume->data[p]/temp1);
				model.data[p][1] *= (_diffractionVolume->data[p]/temp1);
			}
		}

		// Back in real space
		model.fft(FFTW_FORWARD);


		/*
		 *	Maintain running sums
		 */
		if ( (i % _settings->swFreq) > _settings->swFreq/2) {
			strcat(processingString, " ss ");
			for(long p=0; p<nn; p++) {
				shrinkwrapSum.data[p][0] += CPLX_ABS(model.data, p);
				shrinkwrapSum.data[p][1] = 0;
			}
		}

		if ( i > _settings->rsStart) {
			strcat(processingString, " rs ");
			for(long p=0; p<nn; p++) {
				runningSum.data[p][0] += CPLX_ABS(model.data, p);
				runningSum.data[p][1] = 0;
			}
		}

		/*
		 *	Update supoprt
		 */
		if (_settings->sw && i%_settings->swFreq==0 && i>=_settings->swStart && i<_settings->swStop  ) {
			for(long p=0; p<nn; p++) {
				tempModel.data[p][0] = shrinkwrapSum.data[p][0];
				tempModel.data[p][1] = shrinkwrapSum.data[p][0];
				//tempModel.data[p][0] = model.data[p][0];
				//tempModel.data[p][1] = model.data[p][1];
			}
            float swradius = _settings->swStartRadius;
			updateSupport(&tempModel, _settings, swradius, _settings->swThresh);

			// Update support mask
			for(long p=0; p<nn; p++)
				(tempModel.data[p][0]==1)?(support[p]=1):(support[p]=0);

			// Write new support to file
			sprintf(filename,"%s/support%05li.h5",_settings->outdir,i);
			printf("Writing current support to file: %s\n",filename);
			for(long p=0; p<nn; p++)
				tempVolume.data[p] = support[p];
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);

			// Zero running sum
			for(long p=0; p<nn; p++) {
				shrinkwrapSum.data[p][0] = 0;
				shrinkwrapSum.data[p][1] = 0;
			}


		}


		/*
		 *	Save solution every now and then
		 */
		if(i%_settings->savefreq == 0 && i != 0) {
			sprintf(filename,"%s/model%05li-abs.h5",_settings->outdir,i);
			printf("Writing current iterate to file: %s\n",filename);
			for(long p=0; p<nn; p++)
				tempVolume.data[p] = CPLX_ABS(model.data, p);
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);

			sprintf(filename,"%s/model%05li-phase.h5",_settings->outdir,i);
			for(long p=0; p<nn; p++)
				tempVolume.data[p] = CPLX_PHASE(model.data, p);
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);

			sprintf(filename,"%s/model%05li-sum.h5",_settings->outdir,i);
			for(long p=0; p<nn; p++)
				tempVolume.data[p] = CPLX_ABS(runningSum.data, p);
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);

		}


		/*
		 *	Print a line of info for the user
		 */
		end = clock();
		dt = end-start;
		dt /= CLOCKS_PER_SEC;

		time(&endt);
		dt = endt-startt;
		dt /= (i+1);

		printf("%05li: %s --> err=%5.5f R=%5.5f %3.2f sec/iter\n", i, processingString, error, derror, dt);

		FILE	*fp;
		fp = fopen(_settings->logfile,"a");
		if (fp != NULL) {
			fprintf(fp, "%05li: %s --> err=%5.5f R=%5.5f %3.2f sec/iter\n", i, processingString, error, derror, dt);
			fclose(fp);
		}



	}
	// End reconstruction loop


	/*
	 *	Output final solutions
	 */
	sprintf(filename,"%s/model%05li-abs.h5",_settings->outdir,_settings->niterations);
	printf("Writing current iterate to file: %s\n",filename);
	for(long p=0; p<nn; p++)
		tempVolume.data[p] = CPLX_ABS(model.data, p);
	tempVolume.shift(nx/2, nx/2, nx/2);
	tempVolume.writeHDF5(filename);

	sprintf(filename,"%s/model%05li-phase.h5",_settings->outdir,_settings->niterations);
	for(long p=0; p<nn; p++)
		tempVolume.data[p] = CPLX_PHASE(model.data, p);
	tempVolume.shift(nx/2, nx/2, nx/2);
	tempVolume.writeHDF5(filename);


	/*
	 *	Cleanup temp variables
	 */
	printf("End of reconstruction loop\n");

	free(dvMissingData);
	free(support);
}


/*
 *	Update support using Shrinkwrap procedure
 *	Model gets destroyed and replaced with binary mask containing support
 *	(a bit inefficient -> fix later)
 */
void updateSupport(cFFTWdata *model, cSettings *settings, float width, float threshold){

	long	nx = model->nx;

	printf("Updating support with r=%f, thresh=%f\n",width, threshold);
	
	
	/*
	 *	Use temporary array so that model remains unchanged
	 */
	cFFTWdata temp;
	temp.create(nx);
	temp.createFFTWplan(settings->nthreads);

	
	/*
	 *	Take ABS value of current guess
	 */
	for(long p=0; p<temp.nn; p++) {
		temp.data[p][0] = CPLX_ABS(model->data, p);
		temp.data[p][1] = 0.0;		
	}

	// DC value at (0,0,0) always causes problems
	// temp.data[0][0] = temp.data[1][0];
	
	
	
	/*
	 *	Convolve with a Gaussian of FWHM=w
	 */
	float	w, c1;
	float	xx,yy,zz;
	float	rsq, filter;
	long	p;

	w = width/1.664;	 	
	c1 = (-4*PI*w*w/sqrt(2))/(nx*nx);
    //printf("%f\n",c1);
	
	temp.fft(1);
	
	for(long z=0; z<temp.nz; z++){
		(z<temp.nz/2)?(zz=z):(zz=temp.nz-z); 
		for(long y=0; y<temp.ny; y++) {
			(y<temp.ny/2)?(yy=y):(yy=temp.ny-y); 
			for(long x=0; x<temp.nx; x++) {
				(x<temp.nx/2)?(xx=x):(xx=temp.nx-x); 
				
				rsq = (xx*xx+yy*yy+zz*zz);
				filter = exp(c1*rsq);
				
				p = temp.element(x,y,z);
				temp.data[p][0] *= filter;
				temp.data[p][1] *= filter;

			}
		}
	}
	
	temp.fft(-1);

	/*
	 *	Find the maximum value
	 */
	float	gmax = -1e9;
	
	for(long p=0; p<temp.nn; p++) {
		temp.data[p][0] = CPLX_ABS(temp.data, p);
		temp.data[p][1] = 0.0;		

		if(temp.data[p][0] > gmax)
			gmax = temp.data[p][0];
	}


	/*
	 *	Binarise image with everything above threshold = 1
	 */
	float	thresh = threshold*gmax;

	//printf("%g\t%g\t%g\n",gmax, threshold, thresh);
	
	for(long p=0; p<temp.nn; p++) {
		if(temp.data[p][0] >= thresh)
			model->data[p][0] = 1;
		else 
			model->data[p][0] = 0;
		model->data[p][1] = 0;
		
	}
	
	
}



/*
 *	Main loop for doing the iterative phasing
 */
void phaseit(cData3d *diffractionVolume, cData3d *supportVolume, cData3d *guessVolume, cSettings *settings){

	///	A few arrays for keeping track of things
	long	nx = diffractionVolume->nx;
	long	nn = diffractionVolume->nn;
	char	*dvMissingData = (char *) calloc(diffractionVolume->nn, sizeof(char));;
	char	*support = (char *) calloc(diffractionVolume->nn, sizeof(char));
	char	filename[1024];


	if(dvMissingData==NULL || support==NULL) {
		printf("Error allocating temporary arrays\n");
		exit(1);
	}
	
	for(long p=0; p<diffractionVolume->nn; p++) {
		if(diffractionVolume->data[p] >= settings->prNoisefloor || diffractionVolume->data[p] != settings->dvMissingData)
			dvMissingData[p] = 0;
		else
			dvMissingData[p] = 1;
	}

	
	// Where is the support defined?
	for(long p=0; p<nn; p++) {
		if(supportVolume->data[p] == 1)
			support[p] = 1;
		else 
			support[p] = 0;
	}
	
		
	
	
	/*
	 *	Create initial FFT volume with random phases
	 */
	cFFTWdata	model;
	model.create(nx);
	model.createFFTWplan(settings->nthreads);

	float		phase;
	srand(settings->prRandomSeed);
	
	for(long p=0; p<nn; p++) {
		phase = 2*3.141*rand()/RAND_MAX;
		
		if(dvMissingData[p]) {
			model.data[p][0] = 0.0;
			model.data[p][1] = 0.0;		
		}
		else {
			model.data[p][0] = diffractionVolume->data[p]*cos(phase);
			model.data[p][1] = diffractionVolume->data[p]*sin(phase);		
		}
	}

	// Start off in Real space
	model.fft(FFTW_FORWARD);
	
	
	/*
	 *	Memory of last guess (for HIO, etc)
	 */
	cFFTWdata	lastModel;
	lastModel.create(nx);
	lastModel.createFFTWplan(settings->nthreads);
	
	for(long p=0; p<nn; p++) {
		lastModel.data[p][0] = model.data[p][0];
		lastModel.data[p][1] = model.data[p][0];		
	}
	
	
	/*
	 *	Temporary arrays needed for doing stuff later on
	 */
	cData3d		tempVolume;
	cFFTWdata	tempModel;
	cFFTWdata	shrinkwrapSum;
	cFFTWdata	runningSum;
	tempModel.create(nx);
	tempModel.createFFTWplan(settings->nthreads);
	shrinkwrapSum.create(nx);
	runningSum.create(nx);
	tempVolume.create(nx);
	
	for(long p=0; p<nn; p++) {
		shrinkwrapSum.data[p][0] = 0;
		shrinkwrapSum.data[p][1] = 0;		
		runningSum.data[p][0] = 0;
		runningSum.data[p][1] = 0;		
	}

	
	/*
	 *	Actual iteration loop
	 */
	time_t			startt, endt;
	clock_t			start, end;		
	float			dt;

	time(&startt);
	for(long iteration=0; iteration < settings->niterations; iteration++) {
		
		char	processingString[100];

		start = clock();

		/// Start off in real space
		
		/*
		 *	Calculate real space error
		 *  Error is the amplitude outside of support bounds, divided by total amplitude of solution.
		 *	Effectively % of amplitude outside of the support after Pmod but before applying real space constraints
		 *	Calculated by reference to intensity of solution (saves having to take the square root)
		 */
		float	temp1=0;
		float	temp2=0;
		float	error=0;
		for(long p=0; p<nn; p++) {
			temp1 = CPLX_SQR(model.data, p);
			temp2 += temp1;
			if(support[p] == 0)
				error += temp1;
		}
		if(error != 0.0)
			error = error/temp2;
		else 
			error=0.0;

		
		/*
		 *	Periodically do error reduction
		 */
		int	thisMethod = settings->prAlgorithm;
        
		if(settings->prErrorReductionFreq != 0 ){
            if ( (iteration-1)%settings->prErrorReductionFreq == 0 || (iteration-2)%settings->prErrorReductionFreq == 0)
                thisMethod = 1;
        }

		
		
		/*
		 *	Apply real space constraints
		 *	Object must exist within known support constraint in real space
		 *		1 - simple error reduction
		 *		2 - RASR error reduction
		 *		3 - Fienup hypbrid I/O algorithm
		 */
		
		switch(thisMethod) {
			
			// Error reduction
			case 1:		
				strcpy(processingString, "err");
				for(long p=0; p<nn; p++) {
					if(support[p] == 0){
						model.data[p][0] = 0.0;
						model.data[p][1] = 0.0;		
					}
				}
				break;
				
			// RAAR
			case 2:
				strcpy(processingString, "raar");
				for(long p=0; p<nn; p++) {
					if(support[p] == 0){
						model.data[p][0] = (1.0 - 2.0*settings->prBeta)*model.data[p][0] + settings->prBeta*lastModel.data[p][0];
						model.data[p][1] = (1.0 - 2.0*settings->prBeta)*model.data[p][1] + settings->prBeta*lastModel.data[p][1];		
					}
				}
				break;
				
			// HIO
			case 3:
				strcpy(processingString, "hio");
				for(long p=0; p<nn; p++) {
					if(support[p] == 0){
						model.data[p][0] = lastModel.data[p][0] - settings->prBeta*model.data[p][0];
						model.data[p][1] = lastModel.data[p][1] - settings->prBeta*model.data[p][1];
					}
				}
				break;
		}
		
		
		// Remember current real space array for next iteration
		for(long p=0; p<nn; p++) {
			lastModel.data[p][0] = model.data[p][0];
			lastModel.data[p][1] = model.data[p][1];		
		}
		
		

		/* 
		 *	Operations performed in Fourier space
		 */

		// Go to Fourier space
		model.fft(FFTW_BACKWARD);
		for(long p=0; p<nn; p++) {
			model.data[p][0] /= nn;
			model.data[p][1] /= nn;		
		}
		
		
		/*
		 *	Fourier space error metric (R-factor)
		 *	Convention is to calculate R from complex amplitudes (Fobs) 
		 *	Fobs is defined as amplitude of current guess (and not the intensity). 
		 *	Can only be calculated at points where the input data is defined (data bit set).
		 */
		temp1 = 0;
		temp2 = 0;
		float	derror=0;
		float	total=0;
		for(long p=0; p<nn; p++) {
			if(!dvMissingData[p]) {
				temp1 = CPLX_ABS(model.data, p);
				temp2 = diffractionVolume->data[p];
				derror += fabs(temp2-temp1);
				total += temp2;
			}
		}
		derror /= total;	
		
		
		/*
		 *	Apply Fourier space constraint
		 */
		for(long p=0; p<nn; p++) {
			if(!dvMissingData[p]) {
				temp1 = CPLX_ABS(model.data, p);
				if (temp1 < settings->prEpsilon*diffractionVolume->data[p])
					temp1 = settings->prEpsilon*diffractionVolume->data[p];
				model.data[p][0] *= (diffractionVolume->data[p]/temp1);
				model.data[p][1] *= (diffractionVolume->data[p]/temp1);
			}
		}

		// Back in real space
		model.fft(FFTW_FORWARD);

		
		/*
		 *	Maintain running sums
		 */
		if ( (iteration % settings->swFreq) > settings->swFreq/2) {
			strcat(processingString, " ss ");
			for(long p=0; p<nn; p++) {
				shrinkwrapSum.data[p][0] += CPLX_ABS(model.data, p);
				shrinkwrapSum.data[p][1] = 0;
			}
		}
		
		if ( iteration > settings->rsStart) {
			strcat(processingString, " rs ");
			for(long p=0; p<nn; p++) {
				runningSum.data[p][0] += CPLX_ABS(model.data, p);
				runningSum.data[p][1] = 0;
			}
		}
		
		/*
		 *	Update supoprt
		 */
		if (settings->sw && iteration%settings->swFreq==0 && iteration>=settings->swStart && iteration<settings->swStop  ) {
			for(long p=0; p<nn; p++) {
				tempModel.data[p][0] = shrinkwrapSum.data[p][0];
				tempModel.data[p][1] = shrinkwrapSum.data[p][0];		
				//tempModel.data[p][0] = model.data[p][0];
				//tempModel.data[p][1] = model.data[p][1];		
			}
            float swradius = settings->swStartRadius;
			updateSupport(&tempModel, settings, swradius, settings->swThresh);
			
			// Update support mask
			for(long p=0; p<nn; p++) 
				(tempModel.data[p][0]==1)?(support[p]=1):(support[p]=0);
			
			// Write new support to file
			sprintf(filename,"%s/support%05li.h5",settings->outdir,iteration);
			printf("Writing current support to file: %s\n",filename);
			for(long p=0; p<nn; p++) 
				tempVolume.data[p] = support[p];
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);
			
			// Zero running sum
			for(long p=0; p<nn; p++) {
				shrinkwrapSum.data[p][0] = 0;
				shrinkwrapSum.data[p][1] = 0;		
			}

			
		}

		
		/*
		 *	Save solution every now and then
		 */
		if(iteration%settings->savefreq == 0 && iteration != 0) {
			sprintf(filename,"%s/model%05li-abs.h5",settings->outdir,iteration);
			printf("Writing current iterate to file: %s\n",filename);
			for(long p=0; p<nn; p++) 
				tempVolume.data[p] = CPLX_ABS(model.data, p);
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);

			sprintf(filename,"%s/model%05li-phase.h5",settings->outdir,iteration);
			for(long p=0; p<nn; p++) 
				tempVolume.data[p] = CPLX_PHASE(model.data, p);
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);

			sprintf(filename,"%s/model%05li-sum.h5",settings->outdir,iteration);
			for(long p=0; p<nn; p++) 
				tempVolume.data[p] = CPLX_ABS(runningSum.data, p);
			tempVolume.shift(nx/2, nx/2, nx/2);
			tempVolume.writeHDF5(filename);

		}
		
		
		/*
		 *	Print a line of info for the user
		 */
		end = clock();
		dt = end-start;
		dt /= CLOCKS_PER_SEC;

		time(&endt);
		dt = endt-startt;
		dt /= (iteration+1);
		
		printf("%05li: %s --> err=%5.5f R=%5.5f %3.2f sec/iter\n", iteration, processingString, error, derror, dt);

		FILE	*fp;
		fp = fopen(settings->logfile,"a");
		if (fp != NULL) {
			fprintf(fp, "%05li: %s --> err=%5.5f R=%5.5f %3.2f sec/iter\n", iteration, processingString, error, derror, dt);
			fclose(fp);
		}
		
		
		
	}
	// End reconstruction loop
	
	
	/*
	 *	Output final solutions
	 */
	sprintf(filename,"%s/model%05li-abs.h5",settings->outdir,settings->niterations);
	printf("Writing current iterate to file: %s\n",filename);
	for(long p=0; p<nn; p++) 
		tempVolume.data[p] = CPLX_ABS(model.data, p);
	tempVolume.shift(nx/2, nx/2, nx/2);
	tempVolume.writeHDF5(filename);
	
	sprintf(filename,"%s/model%05li-phase.h5",settings->outdir,settings->niterations);
	for(long p=0; p<nn; p++) 
		tempVolume.data[p] = CPLX_PHASE(model.data, p);
	tempVolume.shift(nx/2, nx/2, nx/2);
	tempVolume.writeHDF5(filename);
	
	
	/*
	 *	Cleanup temp variables
	 */
	printf("End of reconstruction loop\n");
	
	free(dvMissingData);
	free(support);
	
}

void Worker::detectMissingData(char* missingData){
	// Where are regions of missing data?
	for(long p=0; p<_diffractionVolume->nn; p++) {
		if(_diffractionVolume->data[p] >= _settings->prNoisefloor || _diffractionVolume->data[p] != _settings->dvMissingData)
			missingData[p] = 0;
		else
			missingData[p] = 1;
	}
}

void Worker::detectSupport(char* support){
	// Where is the support defined?
	for(long p=0; p<_diffractionVolume->nn; p++) {
		if(_supportVolume->data[p] == 1)
			support[p] = 1;
		else
			support[p] = 0;
	}
}

cFFTWdata Worker::createModelWithRandomPhases(char* dvMissingData){
	cFFTWdata	model;
	model.create(_diffractionVolume->nx);
	model.createFFTWplan(_settings->nthreads);

	float		phase;
	srand(_settings->prRandomSeed);

	for(long p=0; p<_diffractionVolume->nn; p++) {
		phase = 2*3.141*rand()/RAND_MAX;

		if(dvMissingData[p]) {
			model.data[p][0] = 0.0;
			model.data[p][1] = 0.0;
		}
		else {
			model.data[p][0] = _diffractionVolume->data[p]*cos(phase);
			model.data[p][1] = _diffractionVolume->data[p]*sin(phase);
		}
	}
	return model;
}
